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1.  Introduction 

Tungsten  heavy  alloys  (WHAs)  are  attractive  candidates  for  use  as  kinetic  energy 
penetrators  (i.e.,  armor-piercing  projectiles)  because  of  their  relatively  large  mass 
density,  high  melting  point,  and  high  strength  at  elevated  rates  of  loading  (Cai  et  al, 
1995).  Previous  experimental  and  numerical  investigations  reported  in  the  literature 
(Weerasooriya  and  Beaulieu,  1993;  Ramesh,  1994,  Subhash  et  ah,  1994a,  b;  Zhou  et 
ah,  1994;  Zhou  and  Clifton,  1997;  Weerasooriya,  1998,  2003;  O’Donnell  and 
Woodward,  2000;  Woodward  and  O’Donnell,  2000)  reveal  a  variety  of  failure  modes 
exhibited  by  WHAs  strained  at  moderate  to  high  rates,  including  brittle  fracture  at 
W-W  grain  boundaries,  matrix-grain  decohesion,  ductile  matrix  rupture,  transgra- 
nular  cleavage  of  W  grains,  and  adiabatic  shear  banding.  Spallation  (material 
separation  due  to  tensile  wave  propagation)  has  been  observed  in  plane  shock  wave 
experiments  on  this  material  (Dandekar  and  Weisgerber,  1999).  Macroscopic 
constitutive  models  (Stevens  and  Batra,  1998;  Wei  et  al.,  2001;  Yadav  et  al.,  2001) 
typically  used  to  represent  the  two-phase  system  in  numerical  simulations  currently 
lack  a  rigorous  description  of  the  kinematics  of  anisotropic  plasticity  (e.g.  plastic 
spin),  the  role  of  crystallographic  orientation  of  constituent  grains,  and  the  effects  of 
crystal  morphology  on  failure  processes.  Microstructural  characteristics  such  as 
grain  shape,  spatial  arrangement  of  grains,  and  local  crystallographic  orientation  are 
known  to  influence  the  response  of  the  material  to  thermomechanical  loading. 
Experimental  (Bruchey  et  al.,  1991,  1992;  Kingman,  1997)  and  numerical 
(Schoenfeld  and  Benson,  1997)  methods  have  demonstrated  a  possible  correlation 
between  lattice  orientation  of  W  single  crystals  and  the  performance  of  such  crystals 
as  kinetic  energy  penetrators.  Wei  et  al.  (2000)  studied  the  influence  of  pre-twisting  of 
polycrystalline  WHA  specimens  on  the  response  under  combined  compressive-shear 
loading  and  found  that  a  certain  degree  of  pre-twisting  of  the  micro  structure 
promotes  adiabatic  shear  localization  at  high  strain  rates.  Such  shear  localization  is 
considered  desirable  in  armor  penetration  events,  since  it  is  thought  to  promote  a 
“self-sharpening”  effect  leading  to  greater  depth  of  penetration  (Magness,  1994; 
Yadav  et  al.,  2001). 

The  particular  material  studied  here  consists  of  relatively  stiff  and  brittle  pure 
tungsten  grains  (BCC),  most  often  nearly  equi-axed  in  shape,  embedded  in  a 
relatively  compliant  and  ductile  matrix  (FCC)  consisting  of  nickel  (50wt.%),  iron 
(25wt.%),  and  tungsten  (25wt.%).  Nominally,  the  composite  micro¬ 
structure  consists  of  90%  pure  W  and  10%  matrix  alloy,  leading  to  an 
overall  weight  distribution  of  93W-5Ni-2Fe.  Fabrication  of  the  composite 
micro  structure  is  conducted  via  isostatic  pressing  and  sintering  of  a  mixture  of  W, 
Ni,  and  Fe  powders,  followed  by  annealing  to  remove  absorbed  hydrogen  and 
then  possible  quenching,  swaging,  and/or  pre-twisting  to  alter  the  dynamic 
mechanical  properties  prior  to  deformation  testing  at  high  strain  rates  (Weerasoor¬ 
iya,  1998;  Wei  et  al.,  2000).  Typical  grain  sizes  are  10-30  pm  for  the  W  crystals  and 
200-500  pm  for  the  matrix  phase  (Zamora  et  al.,  1992;  Zhou,  1993),  meaning 
that  multiple  W  crystals  are  often  embedded  within  each  single  “grain”  of  the 
matrix  phase. 
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Previous  computational  work  by  Zhou  and  co-workers  (Zhou  et  ah,  1994;  Zhou, 
1998a,  b)  elucidated  the  role  of  grain  morphology,  heat  conduction,  strain-  and 
strain-rate  hardening,  and  thermal  softening  on  the  elastoplastic  deformation  and 
shear  localization  of  WHAs  at  nominally  high  applied  rates  of  shear  and  combined 
pressure-shear  loading.  While  grain  and  matrix  phases  were  resolved  discretely  in 
these  calculations,  each  constituent  was  treated  as  an  isotropically-hardening, 
hypothermoelastic-viscoplastic  material,  thereby  neglecting  the  possibly  significant 
effects  of  lattice  rotation  and  anisotropic  strain  hardening  associated  with  the 
evolution  of  crystallographic  texture  within  each  phase.  Thus  the  potential  exists  for 
development  of  more  descriptive  constitutive  models  accounting  for  local  crystal¬ 
lographic  orientation  effects,  models  which  would  simultaneously  engender  more 
accuracy  to  ballistic  simulations  and  support  the  design  of  WHA  microstructures 
tailored  for  enhanced  thermomechanical  properties.  Recent  high  rate  experiments  of 
Weerasooriya  (2003)  indicated  that  tensile  failure  of  WHA  specimens  often  initiates 
via  local  fracture  at  W-W  (i.e.,  grain-grain)  interfaces,  and  less  often  at  interfaces 
between  W  grains  and  the  matrix  phase.  Considering  W-Ni-Fe  alloys  of  various 
compositions  deformed  quasi-statically  and  at  room  temperature.  Woodward  and 
O’Donnell  (2000)  ranked  contributions  of  various  mechanisms  to  total  crack  area,  in 
order  of  percentage  of  potential  area  cracked,  as  follows  in  order  of  descending 
frequency:  W-W  interfacial  failure,  W-matrix  interfacial  separation,  matrix  rupture, 
and  cleavage  fracture  within  W  grains.  However,  in  a  companion  publication, 
O’Donnell  and  Woodward  (2000)  emphasized  that  at  higher  ambient  temperatures 
(>373  K),  the  propensity  for  cleavage  fracture  within  grains  tends  to  increase  at  the 
expense  of  interfacial  fracture.  Zamora  et  al.  (1992)  observed  that  microscopic 
fracture  mechanisms  and  macroscopic  fracture  toughness  measurements  vary  with 
chemical  composition,  grain  size,  and  grain  contiguity.  Failure  strengths  of  interfaces 
within  the  composite  micro  structure  also  depend  strongly  upon  processing 
conditions;  for  example,  slow  cooling  rates  following  post-sintering  annealing  may 
result  in  segregation  of  strength-reducing  impurities  at  interfaces  (German  et  al., 
1984).  In  our  numerical  framework,  we  model  interface  fracture  by  invoking  the 
cohesive  zone  approach  (Barenblatt,  1959;  Dugdale,  1960;  Needleman,  1987;  Xu  and 
Needleman,  1994;  Camacho  and  Ortiz,  1996).  Thus,  unlike  prior  computational 
investigations  that  did  not  consider  fracture  explicitly,  we  are  able  to  capture  a 
failure  mode  that  may  dominate  the  response  of  the  WHA  material  system  when 
subjected  to  a  net  tensile  hydrostatic  stress. 

The  present  work  employs  single  crystalline  plasticity  models  for  each  phase,  as 
discussed  in  Section  2.  Models  for  both  materials  are  embedded  within  a  general 
thermodynamic  framework  applicable  for  describing  any  thermoelastic- viscoplasti- 
cally-deforming  single  crystal.  Finite  deformation,  strain-rate  dependence,  heat 
conduction,  thermal  expansion,  thermal  softening,  and  thermoelastic  coupling  are 
included  in  the  description.  Furthermore,  our  framework  enables  prediction  of  the 
fraction  of  plastic  work  stored  in  the  material  as  residual  elastic  energy  associated 
with  lattice  defects  (i.e.,  dislocations).  Our  treatment  of  this  “stored  energy  of  cold 
work”  can  be  viewed  as  an  extension  to  finite  crystal  plasticity  theory  of  the 
macroscopic,  linearized  elastic-plastic  framework  of  Rosakis  et  al.  (2000),  who 
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demonstrated  the  strain  and  strain-rate  dependence  of  the  fraction  of  plastic  work 
converted  to  stored  energy. 

Cohesive  fracture  modeling  is  discussed  briefly  in  Section  3.  A  temperature-  and 
rate-independent,  linear-softening,  2D  cohesive  finite  element  approach  similar  to 
that  of  Camacho  and  Ortiz  (1996)  is  employed  here,  with  new  fracture  surfaces 
generated  at  interfaces  between  continuum  elements  when  the  traction  on  the 
interface  exceeds  the  intrinsic  fracture  strength  of  the  interface  (a  material 
parameter).  In  our  computational  modeling  scheme,  fracture  surfaces  are  only 
generated  at  interfaces  between  grains  of  the  same  and/or  different  phases;  i.e., 
intragranular  fractures  (cleavage  of  W-grains  and  matrix  rupture)  are  not  included  in 
the  cohesive  zone  description  invoked  in  the  present  work. 

Two-dimensional  finite  element  simulations  of  deforming,  realistic  WHA 
micro  structures  are  conducted,  with  results  given  in  Section  4.  We  investigate 
parametrically  the  effects  of  varying  the  ratio  of  fracture  strengths  of  W-W 
interfaces  and  W-matrix  interfaces,  the  initial  lattice  orientations  of  the  constituent 
W  crystals  and  matrix  phase,  and  the  representation  of  the  micro  structural 
morphology,  the  latter  obtained  from  an  optical  image  of  a  test  sample  of  the 
actual  material.  Each  of  these  aspects  is  shown  to  influence,  to  a  varying  degree,  the 
average  stress-strain  response  and  failure  behavior  of  the  aggregate.  Finally, 
multiscale  issues  important  in  the  construction  of  a  macroscopic  damage  and  failure 
model  for  the  homogenized  material  system  are  mentioned,  following  the  framework 
of  Clayton  and  McDowell  (2003,  2004).  These  issues  include  description  of  the 
contribution  of  damage  entities  (i.e.,  fracture  surfaces)  to  the  average  deformation 
gradient  for  a  representative  volume  of  material  and  the  appropriate  choice  of 
internal  state  variables  capturing  the  evolution  of  damage  and  its  effect  upon  the 
average  stress  supported  by  representative  volume  of  material. 


2.  Crystal  plasticity  formulation 

2.7.  Kinematics 


Let  X  =  x(X,  t)  represent  a  smooth  time-dependent  motion.  A  multiplicative 
decomposition  of  the  deformation  gradient  is  assumed: 


75  If 

f  =  _  —  fef9fP 
5X  “  ’ 


(1) 


where  T,  f®,  and  represent  the  kinematics  of  elasticity  and  rigid-body  rotation, 
thermal  expansion  or  contraction,  and  the  cumulative  contribution  of  moving  crystal 
defects,  respectively.  In  Eq.  (1)  and  henceforward,  juxtaposition  implies  summation 
over  one  set  of  adjacent  indices,  i.e.,  (AB)^^  =  arbitrary  second-rank 

tensors  A  and  B.  In  general,  none  of  T,  f®,  and  F  is  a  compatible  (i.e.,  integrable) 
deformation  gradient  when  considered  individually.  As  shown  in  Fig.  1,  the  elastic 
and  thermal  terms  dictate  the  deformation  of  the  slip  direction  contravariant  vectors 
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Fig.  1.  Deformations,  configurations,  and  slip  vectors. 


and  slip  plane  normal  covariant  vectors 
s(“>  =  rf9s{“t  (2) 

with  the  superposed  “—1”  denoting  the  inverse  operation.  The  velocity  gradient  1 
referred  to  the  current  configuration  is  written  as 

1  =  —  =  ff-'  =  ffe-l  +fepfe-lfe-l  ^ fejejPfp-lfe-lfe-l 

dx  ^ - - - ^  ^ - V - ^ 

=\^  =10  =1P 

where  the  superposed  dot  represents  the  material  time  derivative.  The  thermal 
deformation  is  assumed  isotropic  (Lee  et  al.,  1997),  i.e., 

je  ^  jejo-i  ^  ^4) 

where  6  is  the  temperature  change  measured  from  the  reference  state  and  is  the 
thermal  expansion  coefficient  giving  the  change  in  length  per  unit  current  length  per 
unit  increment  in  6.  The  unit  tensor  is  written  as  1.  The  plastic  velocity  gradient  in 
the  intermediate  configuration  b  of  Fig.  1  is  defined  as  (Rice,  1971;  Asaro,  1983) 

jP  ^  jPjp-l  ^  ^  ^  ^5^ 

a=l 


with  7^“^  the  plastic  shearing  rate  on  slip  system  a,  n  the  number  of  potentially  active 
slip  systems,  and  (g)  the  tensor  product. 
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2.2.  Balance  laws 

Localized  balance  laws  are  prescribed  in  the  current  configuration  as  follows: 
p  +  p  ^r(l)  =  0,  div  (7  +  pb  =  px,  a  =  pe  +  div  q  —  ^r(agl)  =  pr,  (6) 

with  p,  (7,  b,  e,  q,  and  r  the  current  mass  density,  contravariant  Cauchy  stress  tensor, 
body  force  vector  per  unit  mass,  internal  energy  per  unit  mass,  heat  flux  vector  per 
unit  current  area,  and  energy  source  per  unit  mass,  respectively.  Here,  div  denotes 
divergence  with  respect  to  current  coordinates,  e.g.  (i/rq  =  =  daq^  + 

the  Christoffel  symbols  of  the  linear  connection  in  the  spatial  coordinate  system 
satisfying  Igadltc  —  ^cQba  +  ^bQca  ~  ^aQcb^  where  we  use  the  compact  notation  da  = 
d/dx^  and  the  spatial  metric  satisfies  The  trace  operation  for  a  second 

rank  tensor  A  is  denoted  by  tr{A)  =  A\.  The  local  entropy  inequality  is  written  as 
follows,  with  fj  the  time  rate  of  entropy  production  per  unit  mass: 

(7) 

The  Helmholtz  free  energy  per  unit  mass  xj/  is  also  introduced: 

ij/  =  e-et],  (8) 

from  which,  upon  substitution  of  Eqs.  (6)  and  (8)  into  Eq.  (7),  the  entropy  inequality 
becomes 

rr((Tgi)  -  ^  pi'P  +  ^n),  (9) 

with  Vx  the  covariant  derivative  with  respect  to  x  and  •  a  scalar  product  operation 
for  vectors,  i.e.,  a-b  =  dgyV . 


2.3.  Thermodynamic  framework 

We  assume  a  Helmholtz  free  energy  potential  of  the  form 

lA  =  lA(e^0,O,  (10) 

where  the  intermediate  configuration  elastic  strain  2{e^)^^  =  f.a9abft  -  SalS’  with  0^^ 
a  metric  tensor  on  b,  which  in  practice  is  chosen  as  Kronecker’s  delta  for 
simplicity  (Simo  and  Ortiz,  1985),  a  typical  assumption  in  finite  elastoplasticity 
(Clayton  et  ah,  2004).  The  symbol  ^  denotes  a  dimensionless  scalar  internal  variable 
representing  stored  micro-elastic  energy  associated  with  crystal  defects  that  may 
impede  shearing  on  each  slip  system.  Expanding  the  time  rate  of  xj/  and  inserting  the 
result  into  inequality  (9)  yields 

/r(<Tg(P  + 1»  +  IP))  -  +  On),  (1 1) 
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where  the  scalar  product  of  two  second  rank  objects  is  defined  by  A  :  B  =  By,  and 
where  the  subscript  following  the  d-operator  denotes  partial  differentiation  with 
respect  to  the  subscripted  variable.  Additional  algebraic  manipulations  give 


-  pSe'iA)  :  e'  +  (?K®g)ar  -  p(5eiA)  -  pn)Q  - 


where  p  is  the  mass  density  in  configuration  b  and  the  elastic  second  Piola-Kirchhoff 
stress  with  f  =  detT^/det  g/  det  g  =  p/ p  (Teodosiu,  1967). 

The  resolved  Cauchy  stress  for  slip  system  a  is  defined  by  =:  (y(gs^“^ 
Assuming  that  the  independent  variables  in  the  potential  (10)  can  be 
varied  individually  (Coleman  and  Noll,  1963;  Scheidler  and  Wright,  2001),  we  then 
have 


=  pd^exjj, 

(13) 

(Xy 

t]  =  —  tr{ag)  -8g\j/, 

(14) 

y 

=1 

d—l  ^  / 


(15) 


Notice  that  the  first  term  on  the  right  of  Eq.  (14),  denoted  by  /,  arises  as  a 
consequence  of  the  explicit  inclusion  of  thermal  expansion/contraction  in  the 
kinematic  description  (1).  Rearranging  the  energy  balance  in  Eq.  (6)  by  appealing  to 
Eqs.  (13)  and  (14)  leads  to 

n 

p6x  —  pOdeij/  =  pr  —  divq  +  ^  —  p(d^\l/)t.  (16) 

oc—l 


By  defining  the  specific  heat  capacity  c  as 

C  =  8ge  =  d,,eidgt])  =  -0(8ggil/)  +  8^e8gx,  (17) 

where  we  have  invoked  the  differential  of  relation  (8),  and  by  assuming  isotropic 
heat  conduction  in  the  current  configuration  dictated  by  Fourier’s  law  (Zhou  et  al., 
1994): 


q  =  -kv^e, 


(18) 
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with  k(x)  the  thermal  conductivity,  we  can  rewrite  the  localized  energy  balance  (16) 
as 


temperature 

change 


y]  -  pUd^xj/)  -  0id0^\l/))t  +  pddoe'^il/  :  e® 

'  ^  V, _  _ y  V  V 

a—l  ^ 

V  >  energy  of  thermoelastic  coupling 

lattice  defects 

plastic 

dissipation 

+  diviJy^^O)  +  pr  . 

heat  heat 

conduction  supply 


(19) 


2.4.  Constitutive  models 

A  free  energy  potential  per  unit  intermediate  configuration  volume  is  specified  as 

pik  =  ^-e^:  c(e,  0  :  ^  ^(9,  Oe  +  y(0),  (20) 

where  c  and  p  are  the  elastic  modulus  tensor  in  configuration  b  and  an  effective  shear 
modulus,  respectively,  and  /c  is  a  dimensionless,  material-dependent  scalar  parameter 
that  we  assume  is  independent  of  strain  rate  and  temperature.  The  function  y(6)  = 
—c6\n(6/9o)  accounts  for  the  purely  thermal  energy  (Rosakis  et  ah,  2000),  with  Oq  a 
reference  temperature  at  which  y  =  0.  We  shall  later  demonstrate  the  explicit 
relationship  between  k  and  the  fraction  of  plastic  dissipation  converted  to  heat 
energy  (Eq.  (36))  and  will  use  this  relationship  to  select  appropriate  values  of  k. 
From  Eq.  (13)  and  partial  differentiation  of  Eq.  (20),  we  see  that  the  stress  satisfies 
the  linear-hyperelastic  relationship  s®  =  c  :  e®,  or  (s^y^  =  an  adequate 

description  for  the  relatively  small  elastic  strains  that  emerge  in  the  simulations 
conducted  in  the  present  study.  Note  that,  should  the  model  be  used  in  future  work 
for  simulation  of  impact  events  that  may  engender  shockwaves  and  large  volumetric 
strains  in  the  material,  we  suggest  an  augmentation  of  our  description,  for  example 
with  an  equation-of-state  type  formulation  (cf.  Gruneisen,  1926),  to  address  possibly 
nonlinear  relationships  between  pressure,  mass  density,  and  temperature. 

A  power-law  viscoplastic  flow  rule  (Hutchinson,  1976)  is  chosen  to  model  the  time 
rate  of  plastic  deformation  in  each  phase: 

Here  7o  and  m  are  material  parameters,  is  the  slip  resistance  due  to  dislocation 
barriers,  is  the  projected  shear  stress  pulled  back  to  the  intermediate 

configuration  b,  and  sgn(x)  =  x/|x|,  with  sgn(O)  =  1.  Thermal  softening  attributed  to 
increased  dislocation  mobility  at  high  temperatures  is  incorporated  via  the  power- 
law  form  (Klopp  et  al.,  1985) 

=  gt\9/eoy, 


(22) 
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with  the  flow  resistance  at  reference  temperature  Og  (a  material  parameter)  and  p 
a  dimensionless  constant  for  each  material.  We  postulate  the  following  relationship 
between  the  “average”  hardening  on  each  system  at  flxed  reference  temperature  and 
the  internal  variable 

■  ^  “  df)  =  0)  (23) 

a-1 

with  an  initial  yield  stress,  b  the  magnitude  of  the  Burgers  vector,  and  pj  the 
total  dislocation  line  length  per  unit  intermediate  conflguration  volume  associated 
with  shearing  impedance.  The  square-root  dependence  of  flow  stress  on  dislocation 
density  is  a  common  assumption  in  the  plasticity  and  materials  science  literature 
(Taylor,  1934;  Kuhlmann-Wilsdorf,  1985,  1989;  Zikry  and  Kao,  1996;  Kameda  and 
Zikry,  1998;  Ashmawi  and  Zikry,  2003),  as  is  the  assumption  of  linear  dependence  of 
stored  lattice  energy  on  dislocation  density  (Bammann,  2001;  Regueiro  et  al.,  2002; 
Svendsen,  2002)  implied  jointly  by  Eqs.  (20)  and  (23).  The  scalar  proportionality 
factor  a  accounts  for  dislocation  interactions  (Kobytev  et  al.,  1984;  Ashmawi  and 
Zikry,  2003).  Both  lattice  friction  stress  and  effects  of  the  initial  dislocation  density 
are  incorporated  in  g^y  \  the  former  deemed  important  in  characterizing  the  flow 
stress  of  BCC  metals  such  as  tungsten  (Qiu  et  al.,  2001),  the  latter  implying  that  pj  is 
a  measure  of  the  change  in  total  dislocation  density  relative  to  the  initial  state. 

We  now  focus  attention  upon  pure  W  crystals,  for  which  we  permit  slip  in  the 
(111)  close-packed  directions  on  any  of  the  {110}  and  {112}  families  of  planes, 
resulting  in  the  number  of  potentially  active  slip  systems  n  =  24  (Subhash  et  al., 
1994a).  Possible  slip  on  {1  2  3}  planes,  typically  inactive  at  room  temperature  (Argon 
and  Maloof,  1966;  Subhash  et  al.,  1994a),  is  not  represented  in  our  model.  Because 
single  crystalline  W  is  very  nearly  elastically  isotropic,  we  have,  with  the 
dependencies  of  elastic  moduli  on  internal  structure  ^  suppressed  henceforth  for 
simplicity: 

co^uxs  ^  (24) 

with  Lame’s  constant  2  and  contravariant  components  of  the  metric  tensor  on 
conflguration  b.  Evolution  of  slip  resistance  at  reference  temperature  Og  is  specifled 
by  a  hardening-minus-dynamic-recovery  relation  (Armstrong  and  Erederick,  1966; 
Horstemeyer  et  al.,  1999): 

4*)  =  aY)  -  Bg^^^  Y)  (25) 

with  the  interaction  matrix  satisfying 

ql  =  5^p  +  q(l-dp,  (26) 

where  q  is  the  latent  hardening  ratio. 

We  now  discuss  the  constitutive  model  for  W-Ni-Ee  matrix  material,  for  which 
the  number  of  potentially  active  slip  systems  is  chosen  as  /i  =  12.  Dislocation  glide  is 
assumed  to  be  the  dominant  plastic  deformation  mode  for  this  ECC  metal,  occurring 
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in  (110)  close-packed  directions  on  {1  1  1}  planes.  Elastic  isotropy  is  also  assumed 
for  this  phase  (Zhou  et  ah,  1994),  meaning  that  Eq.  (24)  applies,  though  with  elastic 
stiffness  constants  for  the  matrix  substantially  lower  in  magnitude  than  those  for  the 
pure  W.  Strain  rate-  and  temperature-dependent  slip  resistances  in  crystals 
comprising  the  matrix  phase  are  also  specified  via  relations  (25)  and  (26),  albeit 
with  different  values  of  A,  B,  and  q  than  those  used  for  the  pure  W  grains. 


2.5.  Numerical  implementation  and  material  parameters 


Two  idealized  sets  of  thermal  boundary  conditions  are  considered  now  for  the 
purposes  of  calibrating  model  parameters:  isothermal  and  adiabatic.  Note,  however, 
that  heat  conduction  is  fully  taken  into  account  in  the  theoretical  framework  and  in 
the  explicit  finite  element  implementation  discussed  in  later  parts  of  the  current 
work.  A  fully  implicit,  hyperelastic-viscoplastic  algorithm  (Cuitino  and  Ortiz,  1992; 
McGinty,  2001)  is  employed  here  to  integrate  the  elastic-plastic  constitutive 
response.  Let  subscripts  t  and  t  ISt  denote  consecutive  computation  times  in  a 
nonlinear  analysis,  i.e.,  start  and  end  times  in  a  particular  iteration  cycle.  The  slip 
rates  for  a  given  time  increment  spanning  times  t  and  t  +  Ar  are  found  implicitly, 
using  values  of  the  resolved  shear  stress  and  hardening  variables  at  the  end  of  the 
cycle: 


yW  =  to 


^t+St 


(27) 


Since  and  depend  upon  the  solution  variables  an  iterative  procedure  is 
used  to  solve  Eq.  (27).  Notice  that  and  depend  upon  6,  through  Eq.  (22) 
and  the  temperature  dependence  of  elastic  moduli,  Eq.  (24).  For  the  isothermal 
problem,  6  =  0  and  r  =  0  by  assumption,  meaning  that  q  of  Eq.  (18)  varies  spatially 
and  temporally  such  that  the  energy  balance  (19)  is  satisfied  throughout  the  analysis. 
We  also  assume  =  1  for  the  isothermal  case. 

For  the  adiabatic  problem,  we  prescribe  =  0,  r  =  0,  and  initial  conditions 
=  1  and  Ot^o  =  0.  The  temperature  rate  for  a  given  time  increment  spanning 
times  t  and  ^  +  Ans  found  explicitly  using  quantities  at  time  t\ 


where 


^(a)^(a) 


(28) 


(29) 


such  that  1  —  is  the  ratio  of  the  time  rate  of  stored  energy  of  cold  working  to  the 
rate  of  plastic  dissipation  (Taylor  and  Quinney,  1934).  Note  that  while  ^  is  often 
assumed  in  practice  to  acquire  a  fixed  value  for  a  particular  material  regardless  of  the 
temperature-deformation  history,  with  a  value  of  0.8  <  1.0  typical  for  engineering 
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metals  (Kallend  and  Huang,  1984;  Aravas  et  al.,  1990;  Zhou  et  al.,  1994),  we  do  not 
resort  to  such  a  constitutive  assumption  here.  We  do  however  make  the  usual 
assertion  that  the  specific  heat  capacity  for  each  material,  c  =  constant,  a  reasonable 
assumption  for  pure  tungsten  over  the  temperature  range  encountered  in  our 
subsequent  calculations  (Gray,  1972).  The  temperature  at  time  ^  +  is  updated 
simply  via 

0,^^,  =  0,  +  b^t,  (30) 

while  the  thermal  deformation  gradient  at  the  end  of  the  step  is  found  as 

=  exp(ar01Ar)f8,  (31) 

with  exp()  the  matrix  exponential  function.  For  a  given  time  increment  in  an 
adiabatic  analysis,  Eqs.  (27)  are  solved  implicitly  using  values  of  6,  Ot+At,  and 
found  via  (28)-(31).  The  thermoelastic  term  in  Eq.  (19),  upon  assumingy^  =  detf®  ~ 
1  +  3(Xt(0  —  9o)  can  be  rewritten  as 

p9dg,.iP  :  e'  =  9f-^dgis^  :  e*)  -  :  e*) 

%  9f-H[(d0Wl^y^  + 

+  (32) 

The  dependencies  of  the  elastic  moduli  of  the  pure  W  grains  on  temperature  were 
obtained  (to  second  order)  from  the  compilation  of  Yih  and  Wang  (1979): 

de^  =  -3.4  +  0.0065rMPa/°C,  deii  =  -10.3  -  0.0041  rMPa/°C,  (33) 

where  T  =  0  —  273  is  the  temperature  in  degrees  C.  Thermal  dependencies  of  elastic 
moduli  of  the  matrix  phase  were  neglected  in  these  computations,  due  simply  to  lack 
of  experimental  data. 

A  material  point  simulator  was  used  to  generate  average  stress-strain  and 
temperature  histories  for  aggregates  of  one  or  more  grains,  in  order  to  calibrate 

constitutive  model  parameters.  In  these  deformation-controlled  simulations, 

velocities  were  prescribed  as  constant  and  material  accelerations  were  taken  as 
zero.  Also,  Taylor’s  (1938)  assumption  was  invoked  for  polycrystalline  simulations, 
meaning  that  f  was  identical  for  each  grain  in  the  aggregate  throughout  the  time 
history  of  deformation.  For  adiabatic  analyses,  each  grain  maintained  its  own 
temperature,  i.e.,  thermal  interactions  between  crystals  were  forbidden.  Average 
Cauchy  stresses  were  calculated  for  the  aggregate  for  each  time  increment  in  the 
analysis  via  standard  volume  averaging  procedures,  assuming  each  crystal  occupied 
an  equal  volume  in  the  reference  configuration.  We  emphasize  that  the  Taylor 
assumption  and  prescription  of  isothermal  or  adiabatic  conditions  are  enforced  here 
only  for  the  purpose  of  calibrating  material  parameters  in  conjunction  with  a 
material  point  simulator  (i.e.,  a  finite  element  mesh  consisting  of  a  single  element).  In 
later  parts  of  the  present  work,  polycrystalline  aggregates  are  modeled  with 
numerous  finite  elements  per  crystal,  inter-  and  intragranular  interactions  are  taken 
into  account,  inertial  forces  participate,  and  heat  conduction  is  permitted 
throughout  the  domain. 
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Fig.  2  shows  isothermal  stress-strain  curves  for  uniaxial  compression  of  300  pure 
W  grains  of  random  initial  lattice  orientations.  Notice  that  here  we  deform  the 
crystal(s)  in  uniaxial  compression  in  the  2-direction,  i.e.  f\  =  \  —  ht,  with  s  the 
prescribed  nominal  strain  rate.  The  macroscopic  logarithmic  (true)  strain  is  then 
defined  by  s  =  —\n(f\),  while  a  =  f  dV,  with  V  the  total  reference 
configuration  volume  of  the  polycrystal,  equally  partitioned  into  300  grains.  Results 
correspond  to  the  strain-hardening  model  outlined  in  Eqs.  (25)  and  (26),  with 
material  parameters  for  single  crystalline  W  shown  in  column  2  of  Table  1.  Values  of 
A,  and  B  in  Eq.  (25)  were  chosen  based  upon  the  quasi-static  data  for 
poly  crystalline  W  given  by  Dummer  et  al.  (1998),  while  a  relatively  large  value  of 
^  =  1.4  was  selected  for  use  in  Eq.  (26)  following  the  discussion  on  latent  hardening 
in  W  single  crystals  by  Horwath  (1994).  The  strain  rate  and  thermal  sensitivity 
parameters,  m  and  p  respectively,  were  obtained  from  the  crystal  plasticity  model  for 
pure  W  of  Lee  et  al.  (1999).  The  mass  density  in  the  reference  configuration  is  written 
as  Pq,  and  values  of  temperature-dependent  material  properties  (e.g.,  elastic  moduli) 
correspond  to  a  reference  temperature  of  300  K.  Strong  strain  rate  sensitivity  of  the 
flow  stress  is  evident  in  Fig.  2. 

Fig.  3  shows  isothermal  shear  stress-strain  behavior  of  the  pure  matrix  phase, 
calculated  from  the  response  of  300  randomly  oriented  grains.  Shearing  was 
conducted  in  the  1-2  plane,  i.e.,  f  \  =  1  +  The  macroscopic  accumulated  shear 
strain  was  then  found  simply  as  y  =  yt,  while  t  =  V~^  f  a^^dV.  Constants  A  and  B 
in  Eq.  (25)  and  initial  conditions  were  chosen  to  fit  the  quasi-static 

(y  =  10“"^)  torsion  data  of  Zhou  (1993),  and  are  shown  in  column  3  of  Table  1,  along 
with  the  other  material  parameters  for  the  matrix  phase.  The  strain  rate  sensitivity 
exponent  was  assumed  to  be  the  same  as  the  pure  W  grains  following  the 
experimental  data  of  Zhou  (1993).  We  also  selected  yo  to  match  the  value  for  the  pure 
W,  and  used  Taylor’s  (1934)  latent  hardening  assumption  typical  for  FCC  metals. 


Fig.  2.  Average  axial  stress  vs.  compressive  true  strain  for  W,  isothermal  conditions,  300  grains,  Taylor 
model. 
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Table  1 

Material  parameters 


Parameter 

Value  (W) 

Value  (Matrix) 

2 

204  GPa 

137  GPa 

li 

161  GPa 

99  GPa 

Po 

19350  kg/m^ 

9200  kg/m^ 

c 

134J/(kgK) 

382J/(kgK) 

7o 

0.001 

0.001 

m 

20 

20 

q 

1.4 

1.0 

A 

630  MPa 

200  MPa 

B 

1.5 

0.4 

500  MPa 

150  MPa 

p 

-1.5 

-1.5 

eo 

300  K 

300  K 

oct 

5.3(10)“Vk 

i.5(io)-Vk 

k 

160W/(mK) 

100W/(mK) 

a 

0.38 

0.73 

b 

0.318nm 

0.364  nm 

K 

1000 

100 

Fig.  3.  Average  shear  stress  vs.  shear  strain  for  matrix,  isothermal  conditions,  300  grains,  Taylor  model. 


q  =  1.0.  Note  also  that  the  parameters  a  and  k — associated  herein  with  stored  energy 
of  cold  working — need  not  be  specified  in  order  to  conduct  an  isothermal 
stress-strain  prediction. 

For  the  purpose  of  conducting  adiabatic  analyses,  additional  calibrations  were 
needed  to  determine  values  of  a  and  /c  for  each  material,  requiring  particular 
assumptions  on  the  nature  of  stored  energy  of  cold  work.  Manipulating  Eq.  (23)  and 
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then  differentiating  with  respect  to  time  give,  respectively, 


r/im  ^  ^ 


K}i  (Xfin 


Additionally, 


pOideiil,)  =  OdeiKpb^  =  K0{dep){bypy,  =  ^  {gf  -  gf). 

«=i 


Substituting  Eqs.  (34)  and  (35)  into  Eq.  (29)  then  yields 


/?=1 


k{p  +  d{dBfi)) 
(apnf 


E<oS“'-9f)E9o‘' 


Y1  '""f 


(a)  .'.(a) 


a=l 


-1 


(34) 


(35) 


(36) 


_a—l  a— I 

The  dislocation  density  variable  can  be  determined  from  Eq.  (34)  once  a  is  known 
i.e., 

^r=ly,p4>-br>)\.  (37) 


Notice  that  Eq.  (37)  is  a  convenient  relationship  between  the  total  dislocation  density 
pj  and  the  hardness  averaged  over  all  n  potential  slip  systems  at  reference 
temperature  Oq.  The  value  of  parameter  a  can  be  selected  based  upon  experimentally 
obtained  measurements  of  the  dislocation  line  density.  Eollowing  such  a  procedure, 
we  selected  a  value  of  a  by  comparing  simulated  isothermal  results  for  the  dislocation 
density  pj  in  [1  0  0]  and  [1  1  2]  W  crystals  with  the  experimental  data  of  Argon  and 
Maloof  (1966)  on  single  crystals  deformed  quasi-statically  under  uniaxial  tension. 
We  chose  our  value  of  a  for  the  matrix  alloy  by  matching  the  dislocation  evolution  in 
randomly  oriented  300  grain  poly  crystals  of  pure  W  and  matrix  materials.  Nix  and 
Gao  (1998)  used  a  value  of  a  =  0.5  for  copper  and  silver  crystals,  close  to  our  values 
of  0.38  and  0.73  for  W  and  matrix  materials,  respectively,  as  listed  in  Table  1.  It 
should  be  noted  that  since  the  experimental  data  we  used  to  justify  our  choice  of  a 
suffers  from  an  acknowledged  lack  of  precision  (Argon  and  Maloof,  1966),  our 
variable  pj  should  be  viewed  in  the  present  study  as  a  qualitative,  yet  physically- 
based,  measure  useful  for  comparing  the  relative  degree  of  strain  hardening 
accumulated  in  regions  of  W  and  matrix  phases. 

If  data  revealing  evolution  of  P  is  available  from  physical  experiments  (cf. 
Hodowany  et  al.,  2000;  Rosakis  et  al.,  2000),  Eq.  (36)  can  be  inverted  to  deduce  an 
appropriate  value  of  the  parameter  k.  Vice-versa,  if  k  is  known,  then  instantaneous 
values  of  P  can  be  predicted  via  Eq.  (36).  It  is  apparent  from  (20)  and  (36)  that  P  ^  I 
as  /c  ^  0,  meaning  no  energy  of  crystal  defects  is  stored  in  the  lattice  when  k  =  0. 
Caution  must  be  used  when  selecting  a  value  of  k  such  that  j^^O  at  all  times 
(approximately,  in  an  adiabatic  analysis)  in  order  to  satisfy  the  reduced  dissipation 
inequality  (15).  Eig.  4  shows  the  average  value  (^  =  J (^dV  for  randomly 
oriented  polycrystals  of  each  phase  deformed  in  uniaxial  compression,  assuming 
fixed  values  of  k  for  each  analysis  (Table  1)  such  that  0.8<j^<1.0  throughout  the 
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Fig.  4.  Average  cold  work  parameter,  adiabatic  compression,  300  grains,  Taylor  model. 


deformation  history.  Subhash  et  al.  (1994a)  estimated  from  infrared  temperature 
measurements  that  P  approached  unity  in  their  high-rate  compression  tests  on  pure 
W  polycrystals,  while  Zhou  and  co-workers  (Zhou  et  ah,  1994;  Zhou,  1998a,  b)  used 
P  =  0.9  for  matrix  and  pure  W  phases  indynamic  finite  element  calculations.  Notice 
that  instantaneous  values  of  depend  strongly  upon  the  form  of  the  rate  equations 
for  the  slip  resistances  the  latter  which,  upon  consideration  of  Eq.  (34),  can  also 
be  viewed  as  evolution  equations  for  the  defect  density.  For  the  pure  W  phase,  we 
notice  that  P  decreases  with  strain  to  a  minimum  value  around  s  =  0.2,  then  increases 
due  to  the  relatively  large  rate  of  dynamic  recovery  dictated  by  the  choice  of  ^  =  1.5 
in  Eq.  (25).  In  contrast,  for  the  more  steadily  hardening  matrix  phase  (B  =  0.4), 
decreases  continuously  with  increasing  strain.  Please  notice  that  we  include  the  slow- 
rate  adiabatic  data  in  Fig.  4  for  illustrative  purposes,  in  order  to  make  our 
presentation  complete,  since  conventionally  (nearly)  adiabatic  conditions  are 
achieved  only  at  relatively  high  rates  of  loading. 

Figs.  5  and  6  depict  the  stress-strain  behavior  and  temperature  histories  of  both 
materials  under  adiabatic  test  conditions,  where  we  have  used  an  initial  temperature 
of  300  K.  The  volume-averaged  temperature  is  written  as  6.  For  the  matrix  phase, 
upon  using  the  data  of  Zhou  (1993)  to  estimate  an  average  thermal  softening  rate 
dgT,  we  selected  the  same  values  for  thermal  softening  parameters  p  and  0o  (Table  1) 
as  used  for  the  W  phase. 

Significant  heating  and  thermal  softening  is  apparent  in  the  W  material,  especially 
at  the  higher  strain  rate  (s  =  3000/s)  where  the  plastic  dissipation  is  relatively  large. 
Again,  curves  for  slow-rate  adiabatic  tests  are  included  for  illustrative  purposes,  with 
thermal  softening  evident  in  the  W  grains  even  at  an  applied  strain  rate  of 
6  =  10“^ /s.  Comparing  Figs.  3  and  5,  the  relatively  strong  rate  of  strain  hardening 
due  to  dislocation  accumulation  in  the  matrix  phase  is  counteracted  by  thermal 
softening  under  adiabatic  conditions.  Notice  that,  for  a  given  strain  rate,  the 
temperature  rise  in  the  W  grains  tends  to  exceed  that  in  the  matrix,  in  part  because  of 
the  much  larger  specific  heat  parameter  (c)  of  the  latter.  The  melting  temperature  of 
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Fig.  5.  Average  axial  stress  vs.  compressive  true  strain,  adiabatic  compression,  300  grains,  Taylor  model. 


Fig.  6.  Average  temperature  vs.  compressive  true  strain,  adiabatic  conditions,  300  grains,  Taylor  model. 


pure  W  is  3410  K  (Boyer  and  Gall,  1985)  while  that  of  the  matrix  is  approximately 
1750  K  (Zhou  and  Clifton,  1997).  Neither  material  attained  its  melting  point  in  our 
simulations. 


3.  Fracture  modeling 

As  discussed  by  Weerasooriya  (2003),  failure  of  WHA  specimens  subjected  to  high 
rates  of  macroscopic  tensile  deformation  often  initiates  via  local  fracture  at  W-W 
(i.e.,  grain-grain)  interfaces,  and  less  often  at  interfaces  between  W  grains  and  the 
matrix  phase.  The  former  interfaces  are  thought  to  be  weaker  than  the  latter,  based 
upon  visual  examination  of  recovered  specimens  which  demonstrate  that  fracture 
surfaces  are  prone  to  initiate  at  W-W  boundaries  and  then  either  propagate  along 
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W-matrix  interfaces,  or  less  often,  along  cleavage  planes  within  grains  of  either  phase 
(Zamora  et  ah,  1992;  Weerasooriya  et  ah,  1994;  O’Donnell  and  Woodward,  2000; 
Woodward  and  O’Donnell,  2000;  Weerasooriya,  2003). 

Consistent  with  these  experimental  observations,  in  our  simulations  cracks  are 
generated  at  interfaces  between  grains  of  the  same  or  different  phases:  W-W  grain 
boundaries  and  W-matrix  boundaries.  Intragranular  cracking  and  fracture  at 
matrix-matrix  grain  boundaries  are  not  captured,  as  these  failure  modes  were 
observed  less  frequently  in  the  aforementioned  high-rate  experiments  (Weerasooriya 
et  al.,  1994;  Weerasooriya,  2003)  and  in  the  room  temperature,  quasi-static 
experiments  of  Woodward  and  O’Donnell  (2000).  The  dynamic  finite  element 
approach  is  employed,  with  new  fracture  surfaces  generated  at  interfaces  between 
continuum  elements  when  the  traction  acting  upon  the  potential  initiation  site 
exceeds  the  intrinsic  failure  strength  of  the  interface  (a  material  parameter).  Hence, 
fracture  at  interfaces  initiates  when  one  of  the  following  local  stress-based  criteria  is 
attained: 

^  =  ^0,  T  =  To,  (38) 

where  ^  and  t  are  the  resolved  normal  traction  and  shear  traction  on  the  interfacial 
surface,  measured  per  unit  reference  configuration  area,  and  and  tq  are  material 
parameters  specifying  the  normal  and  tangential  (i.e.,  mode  I  and  mode  II)  fracture 
strengths  of  the  interface.  In  the  finite  element  implementation,  duplicate  nodes  are 
generated  along  all  potential  fracture  surfaces  during  the  meshing  stage.  Initially 
coincident  nodes  are  then  constrained  to  share  the  same  velocity  and  temperature 
histories  until  either  of  conditions  (38)  is  reached. 

After  initiation  of  fracture,  the  constitutive  response  of  the  degraded  material  (i.e., 
within  the  cohesive  zones)  at  the  interfaces  is  dictated  by  the  following  set  of  coupled 
irreversible  traction-displacement  relationships  (see  e.g.  Camacho  and  Ortiz,  1996): 

5  =  50^1-^^  (loading), 

5  =  5o  ^  (unloading),  (39) 

T  =  To  ^1  -  ^1  -  sgn(d,)  (loading), 

'  =  (40) 

where  dn  and  dt  are  the  normal  and  tangential  crack  opening  displacements,  dn\  and 
are  the  maximal  values  of  dn  and  dt  achieved  during  prior  loading,  and  Sc  is  a 
material  parameter  describing  the  separation  distance  at  which  the  cohesive  interface 
no  longer  supports  traction  (i.e.,  critical  opening  displacement  for  total  failure).  The 
Macaulay  bracket  is  written  (),  satisfying  (x)  =xVx^0  and  (x)  =  0  Vx<0.  As 
shown  in  Fig.  7(a)  for  pure  tension  and  Fig.  7(b)  for  pure  shear,  unloading 
(reloading)  occurs  linearly  to  (from)  the  origin  in  traction-displacement  space.  Also 
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Fig.  7.  Cohesive  traction-displacement  behavior  for  pure  tension  (a)  and  pure  shear  (b). 


prior  to  attainment  of  fracture  criteria  (38),  heat  conduction  across  interfaces  is 
permitted.  Upon  initiation  and  subsequent  normal  separation  beyond  3c,  however, 
we  enforce  the  null  heat  flux  conditions  S/^9-n  =  0,  where  n  is  the  outward  normal 
vector  to  the  newly-created  fracture  surface. 

It  should  be  noted  that  the  rate-  and  temperature-independent  traction-separation 
laws  described  in  relations  (38)-(40)  were  chosen  because  of  their  wide  acceptance  in 
the  literature  and  their  relative  simplicity,  requiring  few  material-dependent  fracture 
properties.  Particular  choices  of  values  of  these  properties  are  discussed  in  the  next 
Section.  More  elaborate,  and  potentially  more  physically  realistic,  cohesive 
constitutive  models  have  been  advanced  in  the  recent  literature  for  a  variety  of 
material  systems,  and  dynamic  fracture  predictions  are  known  to  exhibit  sensitivity 
to  the  particular  choice  of  cohesive  law  (Falk  et  al.,  2001).  However,  as  will  be 
discussed  more  in  the  next  Section,  the  fracture  properties  of  the  local  interfaces  in 
the  WHA  material  system  are  difficult  to  characterize  experimentally  (Woodward 
and  O’Donnell,  2000),  and  even  the  macroscopic  fracture  toughness  of  the  alloy  is 
not  known  with  great  precision,  as  is  evident  by  the  broad  range  of  values  reported  in 
the  literature  (Zamora  et  al.,  1992).  Regarding  relations  (38),  ratios  of  tq/^o  deviating 
from  unity  were  not  explored  in  the  present  study,  nor  were  mixed-mode  initiation 
criteria.  Such  limiting  assumptions  were  deemed  acceptable  for  the  present  set  of 
uniaxial  tensile  simulations,  in  which  mode  I  fractures  dominated,  though  we 
anticipate  extending  our  framework  in  the  future  to  explore  the  influence  of  mixed¬ 
mode  activation  criteria  on  predicted  fracture  patterns.  Exploration  of  temperature- 
and  strain-rate-dependence  of  the  cohesive  laws  (cf.  Costanzo  and  Walton,  2002), 
and  statistical  variations  in  cohesive  strengths  among  interfaces  (Espinosa  and 
Zavattieri,  2003a,  b;  Zhou  and  Molinari,  2004)  is  beyond  the  scope  of  the  present 
effort,  and  beyond  the  support  of  the  limited  available  experimental  data  of  the 
WHA  material  system  under  consideration.  Parametric  investigations  probing  these 
phenomena,  as  well  as  effects  of  choice  of  various  traction-separation  curve  shapes 
(cf.  Bjerke  and  Lambros,  2003)  and  coupling  schemes  between  normal  and  tangential 
traction-separation  relations  (Ortiz  and  Pandolfi,  1999;  Espinosa  and  Zavattieri, 
2003b)  are  reserved  for  future  work. 
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4.  Finite  element  simulations:  two-phase  microstructure 

4.1.  Numerical  procedures 

The  constitutive  models  presented  in  Sections  2  and  3  for  thermo-elastoplasticity 
of  each  of  the  WHA  phases  and  interfacial  fracture  were  implemented  within  the 
EPIC  dynamic  wave  propagation  finite  element  solver  (Johnson  et  ah,  1997,  2001). 
In  this  approach,  the  equations  of  motion,  i.e.  the  second  of  Eqs.  (6),  are  integrated 
directly  and  explicitly  using  the  algorithm  of  Belytschko  et  al.  (1976).  The 
deformation  gradient  is  then  updated  within  each  element  as 

f,+A,  =  exp(lA0f/,  (41) 

where  1  =  ff“'  is  the  velocity  gradient  that  is  assumed  constant  over  the  time  interval 
{t,  t  +  A^).  The  constitutive  update  proceeds  within  each  element  via  the  methodol¬ 
ogy  discussed  in  Section  2.5.  Contributions  to  the  energy  balance  (19)  from  plastic 
dissipation,  lattice  defects,  and  thermoelastic  coupling  are  calculated,  and  then  the 
temperature  field  is  updated  following  the  explicit  integration  procedure  of  Johnson 
(1981).  Possible  contributions  from  cohesive  elements  to  global  mechanical  and 
thermal  force  vectors  are  accounted  for  just  prior  to  the  enforcement  of  boundary 
conditions  and  initiation  of  the  integration  step  for  the  deformation  field  of  the  next 
cycle.  The  stable  time  increment  for  each  cycle  is  chosen  as  a  small  fraction  of  the 
time  required  for  a  longitudinal  elastic  stress  wave  to  traverse  the  smallest  element  in 
the  grid. 

Results  from  two  different  2D  meshes,  shown  in  Figs.  8(b)  and  8(c),  reconstructed 
from  an  optical  micrograph  of  a  sectioned  micro  structure  obtained  from  an 
undeformed  WHA  sample  (Fig.  8(a)),  will  be  discussed  in  the  present  work.  The 
meshes  consist  of  constant  strain  triangular  elements,  generated  with  the  PPM200F 
software  package  (Langer  et  al.,  2003),  readily  enabling  refinement  along  material 
interfaces  (i.e.,  mesh  refinement  along  potential  cohesive  fracture  surfaces).  The  grid 
in  Fig.  8(c)  was  generated  by  rotating  that  in  Fig.  8(b)  by  90  °,  with  the  latter  a  direct 
reproduction  of  the  image  shown  in  Fig.  8(a).  The  square  domain  is  of  size 
L  =  150  pm,  and  the  area  fraction  of  W-Ni-Fe  matrix  phase  in  each  mesh  is  12.6%. 


(a) 

Fig.  8.  Optical  micrograph  and  finite  element  meshes  of  two-phase  microstructure:  (a)  image,  (b)  mesh  1, 
(c)  mesh  2. 
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Fifty  distinct  W  grains  are  resolved,  with  grids  consisting  of  a  total  of  16664  standard 
3-node  triangular  elements  and  up  to  1656  4-node  cohesive  finite  elements  inserted  at 
potential  fracture  initiation  sites  along  material  interfaces. 

We  enforced  the  following  boundary  and  initial  conditions.  Analyses  were  plane 
strain  tension  in  the  plane,  meaning  that  our  simulations  can  be  thought  to 

represent  columnar  poly  crystals  extended  infinitely  in  the  out-of-plane  direction. 
Please  note  that  out-of-plane  elastic  and  plastic  deformations  were  permitted  (i.e., 
the  3D  material  models  discussed  in  Section  2  were  employed,  with  the  full  number 
of  slip  systems  enabling  out-of-plane  lattice  rotations),  so  long  as  the  total 
deformation  field  remained  planar.  Let  the  lower  and  upper  edges  of  the  domain  be 
denoted  by  =  0  and  =  L,  respectively.  And  let  the  left  and  right  edges  be 
denoted  by  =  0  and  =  L,  respectively.  The  velocity  boundary  conditions  are 
summarized  as 

=  0  along  X^  =  0, 
x'^  =  1.5  m/s  along  X^  =  L, 

=  0  at  (X^  =  L/2,  =  0).  (42) 

Additionally,  force  boundary  conditions  were  applied  along  the  lateral  edges  X^  =  0 
and  X^  =  L  such  that  these  edges  were  constrained  to  remain  straight  and  parallel, 
yet  free  to  contract  upon  extension  of  the  mesh  in  the  -direction.  These  lateral 
edges  were  free  of  shear  traction;  i.e.,  =  0  along  X^  =  0  and  X^  =  L.  The  sides 

were  forced  to  remain  straight  in  order  to  prevent  necking  or  relative  shearing  of 
large  sections  of  the  mesh  due  to  potentially  highly  localized  deformation.  Each 
lateral  edge  was  constrained  in  practice  by  assigning  a  uniform  acceleration  x^  to  all 
nodes  comprising  the  edge,  with  this  value  of  x^  calculated  by  dividing  the  total 
reaction  force  along  the  edge  by  the  total  mass  of  the  nodes  comprising  the  edge,  a 
methodology  similar  to  that  employed  by  Zhou  et  al.  (1994)  for  imposing  periodic 
shear  deformations.  Such  periodicity  was  also  intended  to  instill  the  mesh  with  a 
global  deformation  mode  representative  of  an  actual  element  of  material  within  a 
deforming  test  sample,  i.e.,  a  representative  volume  element  (RVE).  However,  it  is 
noted  that  we  do  not  presume  a  priori  that  our  volume  of  two-phase  material 
contains  a  sufficient  number  of  grains  to  satisfy  the  definition  of  an  RVE  in  a  strict 
sense  (Hill,  1963).  In  order  to  minimize  effects  of  tensile  shock  waves  that  would 
arise  under  impulsive  loading,  the  following  initial  conditions  were  applied 
throughout  the  domain: 

x^  =  (1.5A^/L)m/s  at  r  =  0.  (43) 

Null  heat  flux  boundary  conditions  were  enforced  as 

'  =  0 
X^=L 
=  0 
.X^=L, 


Vx^-n  =  0  along  < 


(44) 
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with  n  the  outward  normal  vector  in  the  spatial  frame,  and  an  initial  temperature 
9q  =  300  K  was  prescribed  uniformly  throughout  the  problem  domain.  Conditions 
(42)  result  in  an  applied  stretch  rate  of  s=  lO'^/s.  Simulations  were  generally 
conducted  for  a  total  duration  of  20  ps,  producing  a  nominal  stretch  of  e  =  0.20  at 
their  conclusion.  Total  execution  times  averaged  960  CPU-hours  divided  among  10 
parallel  processors  on  a  Silicon  Graphics  3900  supercomputer.  Please  notice  that 
tension  simulations  were  of  interest  for  the  study  of  (predominantly)  mode  I  fracture 
in  the  WHA  material  system,  as  this  phenomena  has  not  been  probed  numerically  in 
previous  computational  investigations  of  this  material’s  behavior  (cf.  Zhou  et  al., 
1994;  Zhou,  1998a,  b).  Note  also  that  even  though  tension  tests  are  simulated  here, 
the  constitutive  models  for  plasticity  in  each  phase  of  the  material  were  calibrated 
primarily  using  experimental  compression  and  torsion  data,  as  data  for  tension  is  not 
readily  available  at  higher  strains  and  strain  rates  due  to  the  relatively  brittle  nature 
of  the  material  under  such  loading  conditions. 

Results  from  seven  distinct  simulations  are  discussed  in  the  present  work,  listed  as 
case  numbers  1,2, ...  7  in  Table  2.  For  case  1,  fractures  were  prohibited.  Random 
initial  lattice  orientations  for  the  W  grains  and  matrix  phase  were  assigned.  Each 
pure  W  grain  was  assigned  a  different  initial  lattice  orientation  (i.e.,  set  of  Euler 
angles),  while  the  entire  matrix  phase  was  assigned  the  same  initial  orientation,  such 
that  all  50  pure  W  grains  were  effectively  embedded  within  a  single  “grain”  of  the 
matrix  phase.  In  column  6  of  Table  2,  “1”  refers  to  mesh  1  of  Fig.  8(b).  For  case  2, 
all  parameters  were  identical  to  those  of  case  1,  except  that  fracture  was  permitted 
along  W-matrix  interfaces  only,  with  a  prescribed  strength  for  opening  of  cohesive 
surfaces  at  such  interfaces  of  =  2.0  GPa.  For  case  3,  all  parameters  were 

identical  to  those  of  case  1,  except  that  fracture  was  permitted  along  W-W 
boundaries  only,  with  a  prescribed  strength  for  opening  of  these  cohesive  surfaces  of 
^0  =  To  =  2.0  GPa.  For  case  4,  all  parameters  were  identical  to  those  of  case  1,  except 
that  fracture  was  permitted  along  W-W  boundaries  and  W-matrix  interfaces,  with  a 
uniform  initiation  strength  for  opening  of  potential  cohesive  surfaces  at  all  interfaces 
prescribed  by  ^0  =  ^0  =  2.0  GPa.  For  case  5,  all  parameters  were  identical  for  those 
of  case  4,  except  that  the  initial  orientation  of  the  lattice  of  the  matrix  phase  was 
varied,  as  indicated  by  the  value  “2”  in  column  4  of  Table  2  (an  orientation  with  a 


Table  2 

Numerical  simulations 


Case 

^0,  W-W 
(GPa) 

io,  W-matrix 
(GPa) 

Matrix 

orientation 

W 

orientations 

Mesh 

(morphology) 

1 

00 

00 

1 

1 

1 

2 

00 

2.0 

1 

1 

1 

3 

2.0 

00 

1 

1 

1 

4 

2.0 

2.0 

1 

1 

1 

5 

2.0 

2.0 

2 

1 

1 

6 

2.0 

2.0 

1 

2 

1 

7 

2.0 

2.0 

1 

1 

2 
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significantly  different  Schmid  factor  relative  to  the  other  cases  was  used).  For  case  6, 
all  parameters  were  identical  to  those  of  case  4,  except  that  an  alternate  set  of 
random  initial  lattice  orientations  was  assigned  to  the  50  W  grains,  as  indicated  by 
the  “2”  in  column  5  of  Table  2.  For  case  7,  all  parameters  were  identical  to  those  of 
case  4,  except  that  mesh  2  of  Fig.  8(c)  was  used.  The  same  initial  lattice  orientations 
relative  to  the  global  specimen  axes  were  used  in  cases  4  and  7. 

We  now  discuss  our  choices  of  material  parameters  entering  the  cohesive  relations 
(38)-(40).  For  simulations  in  which  fractures  were  permitted,  nominal  values  of 

=  To  =  2.0  GPa  were  assigned,  in  order  to  roughly  match  peak  stress  levels 
reported  in  experimental  high-rate  macroscopic  data  of  Weerasooriya  (2003),  as  will 
be  discussed  more  later.  Dandekar  and  Weisgerber  (1999)  estimated  from  plate 
impact  experiments  the  spall  threshold  stress  of  WHA  to  lie  between  1.7  and 
2.0  GPa.  In  the  actual  material,  stochastic  variations  in  strength  and  toughness 
among  interfaces  are  of  course  expected  (Woodward  and  O’Donnell,  2000).  Such 
variations  are  not  addressed  in  the  present  work,  but  could  be  readily  implemented 
within  a  cohesive  finite  element  framework  (Zhou  and  Molinari,  2004).  A  uniform 
value  of  dc  =  1.0  pm  was  chosen  for  the  critical  separation  distance  in  the  cohesive 
laws  of  Eqs.  (39)  and  (40),  a  value  thought  to  best  represent  the  effective  fracture 
properties  of  the  WHA  material  system  at  the  length  scale  resolved  by  the  numerical 
discretization  (i.e.,  at  the  scale  of  individual  grains  and  their  interfaces).  Because  the 
WHA  alloys  are  highly  heterogeneous,  with  fracture  behavior  dictated  by 
micro  structural  features  such  as  interfacial  strengths  and  grain  contiguity  and 
influenced  by  large-scale  yielding  (i.e.,  finite  plastic  zones),  difficulties  arise  in 
obtaining  consistent  measurements  of  the  macroscopic  fracture  toughness,  with 
values  of  25<^/c<234MPaVni  reported  by  Zamora  et  al.  (1992)  and  references 
therein.  Furthermore,  the  microscopic  fracture  toughness  of  local  interfaces  is 
expected  to  be  lower  than  that  of  the  homogenized  material,  as  recovered  uniaxial 
test  specimens  reveal  numerous  microcracks  that  do  not  propagate  to  cause 
macroscopic  rupture  (Woodward  and  O’Donnell,  2000).  Assuming  that  the  material 
behaves  linear-elastically  and  neglecting  plastic  dissipation  which  can  dominate 
energy  release  in  the  cohesive  zone  (Rice  and  Wang,  1989),  Woodward  and 
O’Donnell  (2000)  estimated  a  static  fracture  toughness  of  3.4<^/c^7.6MPa^^  for 
microcracking  along  interfaces  within  a  WHA,  while  Diimmer  et  al.  (1998) 
calculated  values  of  0.66  1.5  MPa^^m  for  intergranular  decohesion  in  pure 

polycrystalline  W.  The  mode  I  fracture  energy  corresponding  to  our  choices  dc  = 
1.0  pm  and  =  2.0  GPa  is  =  (l/2)^o^c  =  1-OkJ,  which,  upon  assuming  for 
illustrative  purposes  only,  an  isotropic  linear-elastic  stress-strain  response  for  the 
composite  with  an  effective  elastic  modulus  of  E  =  3.66(10)^  MPa  and  a  Poisson’s 
ratio  of  V  =  0.29  (Dandekar  and  Weisgerber,  1999),  results  in  a  plane  strain  fracture 
toughness  of  Kjc  =  =  20  MPa>/m,  which  we  feel  is  an  acceptable 

compromise  among  the  aforementioned  ranges  of  values  reported  in  the  literature. 
Note  also  that  our  value  of  dc  is  slightly  larger  than  but  of  the  same  order  of 
magnitude  of  the  characteristic  length  of  a  typical  bulk  finite  element  within  the 
mesh,  which  we  found  of  adequate  resolution  to  facilitate  stable,  converged 
solutions.  Cohesive  elements  were  deemed  of  sufficiently  small  size  to  resolve  the 
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process  zone,  the  length  of  which  may  be  liberally  estimated  again  assuming  a  linear 
elastic  bulk  response  (cf.  Rice,  1968;  Espinosa  and  Zavattieri,  2003b). 

4.2.  Numerical  results 

Contour  plots  of  solution  field  variables  corresponding  to  cases  1-7  are  given  in 
Figs.  9-15,  respectively,  all  corresponding  to  ^  =  10  ps,  or  an  applied  tensile  strain  of 
a  =  0.10.  In  part  (a)  of  each  figure,  the  effective  deviatoric  stress  Oe  referred  to  the 
spatial  frame  is  defined  by 

Ge  =  \/3/2((7  —  l/3rr((7)l)  :  (a  —  l/3^r((y)l),  (45) 

while  the  scalar  effective  plastic  strain  &p  of  part  (b)  of  each  figure  is  found  by 

Sp^  J  v'(2/3)dP  :  d^dt,  (46) 

where  2d^  =  F  +  is  the  rate  of  plastic  stretching  referred  to  the  spatial 
configuration.  Part  (c)  of  each  figure  shows  the  absolute  temperature  9,  part  (d) 
gives  the  dislocation  density  pj  introduced  in  Eq.  (23),  and  part  (e)  illustrates 
instantaneous  values  of  the  plastic  work-to-heat  conversion  parameter  P  of  Eqs.  (29) 
and  (36). 

Consider  first  the  damage-free  case  1,  with  contours  given  in  Fig.  9.  From  Fig. 
9(a),  we  notice  that  the  W  grains  generally  support  higher  effective  stress  than  the 
W-Ni-Fe  matrix  phase.  From  Fig.  9(b),  the  opposite  is  true  regarding  effective 
plastic  strain,  with  the  matrix  phase  accommodating  more  of  the  total  applied 
stretch.  Also,  as  seen  in  Fig.  9(c),  the  temperature  rise  is  greatest  in  regions  of  the 
matrix  phase  relative  to  that  in  the  W  grains,  as  a  result  of  relatively  greater  plastic 
deformation  and  subsequent  conversion  to  heat  energy  in  the  former.  Clearly,  over 
the  short  duration  of  the  simulation  (t  =  10  ps  in  Fig.  10),  heat  energy  does  not  have 
time  to  diffuse  uniformly  throughout  the  material  via  conduction.  As  seen  in  Fig. 
9(d),  dislocation  activity  is  also  concentrated  in  matrix  regions,  as  the  ductile  phase 
undergoes  more  straining,  more  strain  hardening,  and  more  dislocation  accumula¬ 
tion  than  the  W  phase.  Fig.  9(e)  shows  P,  which  can  be  thought  of  as  a  snapshot  of 
dislocation  accumulation  rates  at  the  selected  time  of  t  =  10  ps.  The  greater  the  value 
of  at  a  particular  time  instant,  the  lesser  the  rate  of  local  strain  hardening  and 
associated  dislocation  accumulation,  and  the  greater  the  fraction  of  plastic 
dissipation  converted  into  a  temperature  rise.  Regions  of  relatively  high  and  low 
values  of  are  not  restricted  to  one  particular  phase  in  the  material,  indicating  that 
plastic  straining  and  dislocation  accumulation  are  occurring  non-uniformly 
throughout  the  sample  at  t  =  10  ps. 

We  now  discuss  the  results  for  the  elastic-plastic-fracture  simulations  shown  in 
Figs.  10-15  collectively  and  relative  to  those  of  Fig.  9.  Considering  the  stress 
contours  (part  (a)  of  each  figure),  it  is  clear  that  stresses  are  reduced  on  average  as  a 
result  of  microcracking,  though  local  stress  concentrations  in  excess  of  Ge  =  3.0  GPa 
are  present  in  all  cases.  The  net  stress  reduction  is  greatest  in  cases  4-7,  which  permit 
fracture  at  all  interfaces,  as  opposed  to  cases  2  and  3,  in  which  fracture  is  only 
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Fig.  9.  Contours  of  results  at  applied  stretch  of  a  =  0.10  for  case  1:  (a)  effective  stress  Og  (GPa),  (b) 
effective  plastic  strain  Sp,  (c)  temperature  9  (K),  (d)  dislocation  density  pp  (xlO“^/cm^),  (e)  heat 
dissipation  parameter  p. 


permitted  at  W-matrix  or  W-W  interfaces,  respectively.  Homogenized  stress-strain 
curves  for  the  material  samples  will  be  discussed  more  later.  Comparing  Fig.  9(b) 
with  Figs.  10(b)-15(b),  damage  in  the  form  of  microcracking  accommodates  much  of 
the  inelastic  deformation  that  would  otherwise  be  accommodated  by  plastic  straining 
in  the  matrix  phase.  Contours  of  effective  plastic  strain  in  Figs.  10(b)-15(b)  are 
calibrated  to  show  a  maximum  value  of  Sp  =  2.0,  with  all  regions  exhibiting  strains 
greater  than  this  value  red  in  color.  In  the  calculations  for  cases  2-7,  we  imposed  an 
element-based  failure  criterion  in  addition  to  the  cohesive  fracture  model,  removing 
severely  distorted  elements  from  the  global  mechanical  and  thermal  force  balances 
and  the  stable  time  step  estimation  calculation  upon  attainment  of  very  large  local 
strains,  specifically  when  e^>10.0.  This  criterion  permitted  us  to  continue  the 
simulations  to  large  applied  strain  levels  and  also  allowed  some  cracks  to  propagate 
through  regions  of  bulk  material.  This  element  failure  criterion  is  justified  by  the  fact 
that  the  most  severely  distorted  elements  tend  to  have  undergone  significant  thermal 


J.D.  Clayton  /  J.  Mech.  Phys.  Solids  53  (2005)  261-301 


285 


plattic_sttaiB 

lOOOetOO 

IJOnwflO 

l.OOOe+00 

5.<inne^i 

l.272e-Q5 


IJ93C403 

I.195e4^ 

7.9«e402 

3.9«3e+02 

I.2l4e-I6 


Fig.  10.  Contours  of  results  at  applied  stretch  of  a  =  0.10  for  case  2:  (a)  effective  stress  Og  (GPa),  (b) 
effective  plastic  strain  Sp,  (c)  temperature  6  (K),  (d)  dislocation  density  pp  (xl0“^/cm^  ),  (e)  heat 
dissipation  parameter  p. 


softening  as  a  result  of  intense  plastic  deformation,  and  hence  carry  relatively  little 
stress  prior  to  their  removal.  This  element  failure  criterion  was  not  used  for  case  1. 
As  is  clear  from  parts  (b),  (c),  and  (d)  of  Figs.  10-15,  intense  plastic  strain, 
temperature  rise,  and  dislocation  accumulation  are  concentrated  in  the  matrix  phase 
and  within  relatively  distorted  regions  of  the  W  grains  in  the  vicinity  of  failed 
interfaces.  From  part  (e)  of  Figs.  10-15,  P  approaches  unity  throughout  most  of  the 
domain  in  each  case,  although  scattered  regions  with  relatively  low  values  of  P, 
indicating  instantaneous  strain  hardening  activity,  are  evident  in  each  figure. 

Different  microcrack  patterns  and  fracture  paths  are  evident  for  each  of  cases  2-7, 
as  is  clear  from  Figs.  10-15.  The  most  notable  differences  arise  in  cases  2  and  3  (Figs. 
10  and  11,  respectively)  relative  to  cases  4-7,  Figs.  12-15.  Potential  fracture  paths  are 
impeded  in  cases  2  and  3  relative  to  cases  4-7,  because  of  the  limitation  in  the  former 
of  crack  propagation  to  either  W-matrix  or  W-W  boundaries  exclusively.  In  cases 
4-7,  propagation  of  damage  is  less  inhibited  and  mode  I  crack(s)  are  able  to 
propagate  more  readily  across  the  entire  domain.  Thus  our  simulated  cases  4-7  most 
closely  replicate  the  hypotheses  of  Woodward  and  O’Donnell  (2000)  and 
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Fig.  11.  Contours  of  results  at  applied  stretch  of  a  =  0.10  for  case  3:  (a)  effective  stress  Og  (GPa),  (b) 
effective  plastic  strain  Sp,  (c)  temperature  6  (K),  (d)  dislocation  density  pp  (xlO“^/cm^),  (e)  heat 
dissipation  parameter  p. 


Weerasooriya  (2003),  who  suggested  that  cracks  initiate  primarily  at  the  weakest 
W-W  grain  boundaries  and  grow  subsequently  by  following  trajectories  of  minimum 
resistance,  in  our  simulations  along  interconnected  W-W  and  W-matrix  interfaces. 
Comparing  Fig.  11  with  Fig.  12,  it  is  clear  that  limiting  crack  growth  to  W-W 
interfaces  alone  severely  restricts  macro-crack  propagation  relative  to  the  case  when 
both  W-W  and  W-matrix  fractures  are  permitted  simultaneously.  If  we  view  our 
aggregate  of  material  as  statistically  representative,  this  result  suggests  that  W- 
matrix  interface  failures  (or  intragranular  failures)  must  take  place  for  macro-crack 
propagation  in  the  actual  material,  as  the  W-W  boundaries  alone  cannot  account  for 
propagation  across  more  than  a  few  dozen  micrometers,  as  limited  by  the  occurrence 
of  continuously  linked  W-W  grain  boundaries.  Also  from  Figs.  12-14,  we  notice  that 
random  changes  in  lattice  orientation  of  W  grains  or  the  matrix  phase  do  affect  the 
distribution  of  damage,  stress,  and  deformation  within  each  sample,  though 
dominant  cracks  with  similar  paths  traverse  the  upper  and  lower  thirds  of  each 
mesh  in  all  of  these  cases. 
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Fig.  12.  Contours  of  results  at  applied  stretch  of  s  =  0.10  for  case  4:  (a)  effective  stress  Og  (GPa),  (b) 
effective  plastic  strain  Sp,  (c)  temperature  9  (K),  (d)  dislocation  density  pp  (xlO“^/cm^),  (e)  heat 
dissipation  parameter  p. 


Table  3  provides  the  time  instant,  cohesive  element  number,  and  interface  type 
corresponding  to  the  initiation  of  damage  in  each  simulation.  Data  in  Table  3 
correspond  to  satisfaction  of  initiation  conditions  (38)  and  not  total  stress  relief  in 
the  cohesive  zones.  We  see  that  damage  initiation  occurs  very  early  in  each  case 
(t  ^  0.3  ps,  at  an  applied  stretch  of  a  ~  0.3%),  and  always  at  a  W-W  interface  when 
possible  (recall  that  W-W  interface  fractures  were  prohibited  for  case  2). 
Significantly,  first  fracture  initiated  at  W-W  interfaces  even  when  all  interfaces 
were  given  the  same  strength  of  2.0  GPa,  as  in  cases  4-7.  This  is  expected  considering 
our  prior  observation  that  the  W  grains  tend  to  carry  higher  stresses  than  the  matrix 
material,  meaning  that  stresses  should  be  relaxed  at  typical  W-matrix  interfaces 
relative  to  W-W  interfaces.  Since  fracture  occurred  at  the  same  location  and 
virtually  the  same  time  instant  for  each  of  cases  3-6,  we  deduce  that  the  deformation 
at  this  point  throughout  the  mesh  was  nearly  (thermo)elastic,  since  elastic  isotropy 
was  enforced  for  each  phase.  However,  since  the  first  fracture  in  case  5  took  place 
slightly  later  than  that  of  cases  3,  4,  and  6,  some  minor  plastic  deformation  probably 
occurred  in  the  matrix  material  preceding  the  first  fracture.  Recall  that  case  7  is  a 
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Fig.  13.  Contours  of  results  at  applied  stretch  of  a  =  0.10  for  case  5:  (a)  effective  stress  Oe  (GPa),  (b) 
effective  plastic  strain  Sp,  (c)  temperature  9  (K),  (d)  dislocation  density  pp  (xlO“^/cm^),  (e)  heat 
dissipation  parameter  p. 


different  mesh  than  cases  2-6,  hence  a  different  cohesive  element  failed  first  in  the 
former  simulation,  as  would  be  expected. 

Table  1  and  Figs.  10-15  point  to  design  protocols  for  W  alloys  that  could  result  in 
improved  resistance  to  failure.  It  is  suggested  that  by  improving  the  fracture 
strengths  of  the  W-W  interfaces,  or  by  reducing  their  relative  frequency,  fracture 
may  be  delayed.  Such  a  conjecture  agrees  with  the  experiments  of  Weerasooriya  et  al. 
(1994),  who  were  able  to  achieve  increases  in  strain-to-failure  by  reducing  the 
contiguity  of  W  grains  (and  thereby  the  frequency  of  W-W  boundaries)  in  WHAs 
deformed  in  high  rate  torsion  tests  in  a  Kolsky  bar  apparatus.  Furthermore,  our 
results  suggest  that  by  increasing  the  strength  of  one  type  of  interface  relative  to  the 
other,  and  thereby  limiting  fracture  to  just  a  single  mechanism  (i.e.,  only  one  type  of 
interface,  and  no  intragranular  failure),  macro-crack  propagation  may  be  severely 
impeded. 

Fig.  16  presents  average  effective  stress-strain  data  for  each  of  the  simulations, 
de  =  f  (JedV,  with  V  the  total  volume  of  the  problem  domain  in  the  reference 
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Fig.  14.  Contours  of  results  at  applied  stretch  of  a  =  0.10  for  case  6:  (a)  effective  stress  Og  (GPa),  (b) 
effective  plastic  strain  Sp,  (c)  temperature  6  (K),  (d)  dislocation  density  pp  (xl0“^/cm^),  (e)  heat 
dissipation  parameter  p. 


configuration.  In  the  legend  of  Fig.  16,  case  labels  of  each  simulation  consistent  with 
the  identifiers  listed  in  Table  2  are  indicated  in  parentheses.  Data  for  the  damage-free 
case  1  is  given  only  over  the  range  0.0 0.10,  since  an  element  failure  criterion  was 
not  used  in  that  particular  simulation.  The  simulation  was  halted  at  this  strain  level 
due  to  highly  distorted  elements  in  the  vicinity  of  a  localized  deformation  zone 
sustaining  significant  heating.  For  the  remaining  cases  2-7,  data  is  shown  over  the 
applied  strain  range  0.0^a<0.20.  For  case  1,  a  peak  stress  of  de  =  2.55 GPa  occurs 
at  an  applied  stretch  of  a  =  0.024,  beyond  which  the  average  stress  decreases  as  a 
result  of  temperature  rise  and  commensurate  thermal  softening.  Average  stress 
values  are  significantly  reduced  in  cases  2-7  relative  to  case  1  as  a  result  of  damage 
(i.e.,  material  separation  at  cohesive  interfaces).  For  case  2,  in  which  fracture  is 
restricted  to  W-matrix  interfaces,  a  peak  average  stress  of  de  =  1.98  GPa  is  attained 
at  an  applied  stretch  of  e  =  0.012.  For  case  3,  with  fracture  allowed  at  W-W 
interfaces  only,  a  peak  stress  of  de  =  1.93  GPa  is  reached  at  a  =  0.010.  Cases  4-7,  in 
which  all  interfaces  support  potential  fracture,  exhibit  very  similar  average 
stress-strain  behavior,  with  maximum  stresses  of  de  ~  1.52  GPa  attained  at  a  stretch 
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Fig.  15.  Contours  of  results  at  applied  stretch  of  s  =  0.10  for  case  7:  (a)  effective  stress  Og  (GPa),  (b) 
effective  plastic  strain  Sp,  (c)  temperature  9  (K),  (d)  dislocation  density  pp  (xl0“^/cm^),  (e)  heat 
dissipation  parameter  p. 


Table  3 

Times  and  locations  of  initial  fractures 


Case 

Time  (ps) 

Element  # 

Interface  type 

2 

0.3375 

146 

W-matrix 

3 

0.3122 

1309 

W-W 

4 

0.3122 

1309 

w-w 

5 

0.3132 

1309 

W-W 

6 

0.3122 

1309 

w-w 

7 

0.3039 

1227 

w-w 

of  a  ^  0.008  in  each  simulation.  From  these  results  we  arrive  at  the  rather  obvious 
conclusion  that  increasing  the  available  fracture  sites  in  the  sample  results  in  a 
greater  reduction  in  effective  stiffness  and  the  attainment  of  a  peak  average  stress  at 
a  lower  value  of  applied  strain.  Also  notice  from  the  similarity  of  curves  for  cases  4-7 
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Fig.  16.  Average  effective  stress,  Og,  vs.  applied  axial  stretch,  s. 


that  grain  orientation  and  morphology  only  mildly  influence  the  homogenized 
stress-strain  behavior,  lending  us  confldence  that  our  sample  of  material  is  to  a 
certain  degree  statistically  representative  with  regards  to  effective  mechanical 
stiffness,  though  more  tests  should  be  conducted  to  fully  justify  such  an  assertion. 
Weerasooriya  (2003)  performed  quasi-static  and  higher  rate  Kolsky  bar  tension  tests 
on  WHA  samples  and  measured  the  effective  stress-strain  response.  At  the  highest 
rate  conducted,  a  =  750/s,  Weerasooriya  (2003)  predicted  a  peak  uniaxial  stress  of 
1.65  GPa  (corresponding  to  de  =  1.35  GPa)  at  a  strain  of  a  ~  0.01,  followed  by 
failure  of  the  specimen  by  necking  (i.e.,  tensile  instability).  Furthermore,  peak 
stresses  dropped  and  applied  strains  at  peak  stress  increased  as  the  applied  strain  rate 
was  decreased  in  that  series  of  experiments.  It  is  difficult  to  compare  the 
homogenized  stress-strain  behavior  from  our  simulations  to  the  macroscopic  curves 
of  those  experiments  (Weerasooriya,  2003)  directly,  since  different  sample  sizes, 
strain  rates,  and  boundary  conditions  were  enforced  in  the  experiments  and 
simulations.  However,  the  similarity  in  the  range  of  peak  stress  values  attained  in 
experiment  and  simulation  provide  confidence  in  our  selection  of  material 
parameters  used  to  model  plasticity  and  fracture.  Notice  also  that  the  average 
effective  stresses  in  the  fracture  simulations  do  not  relax  completely  to  a  value  of  zero 
upon  attainment  of  a  =  0.20.  Values  of  the  average  stress  de  are  influenced  by  local 
stress  concentrations  near  damage  entities,  residual  stresses  due  to  intergranular 
incompatibility,  and  out-of-plane  and  transverse  stresses  due  to  the  imposed 
macroscopic  boundary  conditions.  Fluctuations  in  stresses  are  also  induced  by 
inertial  effects,  though  these  are  smoothed  substantially  by  the  averaging  process. 
However,  as  will  be  shown  later  in  Section  4.3,  the  average  tensile  stresses  in  the 
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direction  of  elongation  do  relax  towards  zero  magnitude  as  a  result  of  microcracking 
transverse  to  the  loading  direction. 

Average  temperatures  6  for  the  problem  domain  (i.e.,  for  the  entire  FE  mesh)  are 
shown  in  Fig.  17,  with  9  =  f  OdV.  The  observed  trend  is  a  reduction  in  9  with 
increasing  occurrence  of  damage  at  interfaces,  upon  comparing  values  of  9  for  cases 
1-3  with  those  for  cases  4-7  at  particular  time  instances  throughout  the  deformation 
history.  This  phenomenon  is  easily  explained:  microcracks  in  cases  4-7  accom¬ 
modate  much  of  the  deformation  accommodated  by  plastic  strain  in  cases  1-3, 
leading  to  less  cumulative  plastic  dissipation  converted  to  heat  in  the  former. 
Interestingly,  for  small  magnitudes  of  applied  stretch  limited  to  a  <0.05,  rates  of 
average  temperature  increase  8^9  in  cases  2  and  3  exceed  the  rate  of  increase  in  the 
damage  free  case  1,  as  a  result  of  extreme  localized  plastic  deformation  in  the  vicinity 
of  damaged  zones  in  cases  2  and  3.  Such  regions  of  extreme  temperature  rise  (locally, 
0>1OOOK)  are  not  observed  in  cases  4-7,  as  cracks  can  propagate  more  freely 
because  fracture  is  active  at  both  types  of  grain  boundary  interfaces,  and  fewer 
ligaments  of  localized  plastic  deformation  and  intense  heating  are  formed. 

Fig^  18  depicts  the  average  fraction  of  plastic  dissipation  converted  to  temperature 
rise,  P  f  jSdV.  Note  that  0.9<j^<1.0  throughout  the  duration  of  each 

simulation.  Cases  1-3  achieve  lower  values  of  P  throughout  the  deformation  history 
than  cases  4-7,  since  fracture  suppresses  plasticity  and  dislocation  accumulation  in 
the  latter,  especially  in  the  W  grains.  Recall  that  the  W  grains  comprise  most  of  the 
volume  of  the  sample  (87.4%)  and  thus  have  the  strongest  influence  on  volume- 
averaged  quantities  such  as  p.  Rates  of  decrease  —diP  for  cases  2  and  3  exceed  those 
for  case  1  at  low  strain  levels,  again  as  a  result  of  the  emergence  of  ligaments  of 
highly  strained  material  in  the  vicinity  of  damaged  zones  where  cracks  are  unable  to 


-  N o  fracture  ( 1 ) 

.  W-matrix  fracture  only  (2) 

- W-W  fracture  only  (3) 

- W-W  and  W-matrix  fracture  (4) 

- W-W  and  W-matrix  fracture  (5) 

- W-W  and  W-matrix  fracture  (6) 
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Fig.  17.  Average  temperature,  6,  vs.  applied  axial  stretch,  s. 
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- W-W  and  W-matrix  fracture  (6) 

- W-W  and  W-matrix  fracture  (7) 


Fig.  18.  Average  dissipation  fraction,  p,  vs.  applied  axial  stretch,  s. 


propagate  freely  due  to  infinite  cohesive  strengths  of  one  type  of  grain  boundary 
interface. 


4.3.  Macroscopic  damage  modeling 


An  important  objective  of  the  present  research  effort  is  support  of  the  construction 
of  a  physically  realistic  macroscopic  model  of  the  kinematics,  thermodynamics,  and 
kinetics  of  fully  anisotropic  damage  for  subsequent  numerical  implementation. 
Briefly  addressing  the  kinematics  of  such  a  model,  the  net  deformation  gradient  F  for 
a  damaged  material  element  containing  k  internal  surfaces  can  be  decomposed  as 
(see  derivation  based  on  the  generalized  Gauss’s  theorem  in  Clayton  and  McDowell, 
2003) 


F  =  4  /x(g)Nd5'  =  ^  /  idV  +  l-S"  f  x®  O =  F -F F**, 
V  Js  y  Jv  V  Jsm 


(47) 


external 

boundary 


average  material 
deformation,  F 


internal  surfaces 
(damage),  F*^ 


where  S  is  the  external  surface  of  the  volume  element  of  referential  volume  F,  with 
unit  surface  normal  vector  N  (referred  to  the  reference  configuration)  and  surface 
spatial  coordinates  x.  In  Eq.  (47),  F  is  the  volume-averaged  deformation  gradient 
contribution  from  the  intact  material,  and  F*^  is  the  contribution  from  internal 
surface  discontinuities  associated  with  damage  entities  (e.g.  cracks,  voids,  and  shear 
bands),  each  with  referential  surface  reference  normal  vector  and  current 
coordinates  The  premise  of  Eq.  (47)  for  a  polycrystal  undergoing  intergranular 
fracture  is  illustrated  in  Fig.  19. 
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Fig.  19.  Deformation  and  damage  in  a  polycrystalline  volume  element. 


-  W- matrix  fracture  only  (2) 

■  ■■  ■  W-W  fracture  only  (3) 

- W-W  and  W-matrix  fracture  (4) 

—  —  —  -  W-W  and  W-matrix  fracture  (5) 

- W-W  and  W-matrix  fracture  (6) 

- W-W  and  W-matrix  fracture  (7) 


Fig.  20.  Damage  accommodation,  A,  vs.  applied  axial  stretch,  s. 


Shown  in  Fig.  20  is  a  measure  of  the  net  accommodation  of  damage  in  the  tensile 
direction,  for  the  material  volume  elements  (i.e.,  aggregates  of  W  grains,  matrix 
material,  and  cohesive  fracture  surfaces)  in  our  numerical  simulations  2-7,  defined  as 

p'dl 

(48) 

Notice  that  at  applied  strains  a  >0.06,  damage  accommodates  more  than  90%  (i.e., 
yl>0.90)  of  the  total  tensile  stretch  in  cases  4-7,  simulations  in  which  cracks  were 
permitted  to  open  at  W-W  and  W-matrix  boundaries.  In  cases  2  and  3,  in  which 
fractures  were  restricted  to  W-matrix  or  W-W  interfaces  alone,  respectively,  damage 
accommodation  is  limited  to  a  maximum  of  yl<0.80  at  applied  strains  greater  than 
a  =  0.15,  still  a  considerably  large  fraction  of  the  tensile  component  of  the  total 
prescribed  deformation  gradient  F. 
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Accompanying  Eq.  (47)  is  the  definition  of  the  net  nominal  stress  tensor  S,  work 
conjugate  to  the  time  rate  of  the  net  deformation  gradient,  F  (Clayton  and 
McDowell,  2004): 


sdF- 


traction  carried  by 
external  boundary 


stress  in  _ 
material,  S 


-E 


f  : 


X®  0 


traction  carried  by  internal  surfaces,  S** 


+  ^  [  PoX0(x-b)dF  =  S  +  S‘'  +  S‘, 
E  Jv 


local  inertia  and  body  forces,  S* 


(49) 


where  s  =  is  the  local  nominal  stress  in  the  material  (transpose  of  the  first 
Piola-Kirchhoff  stress),  t  is  the  traction  vector  per  unit  reference  area  along  S,  and 
and  are  reference  surface  coordinates  and  traction  vector,  respectively, 
supported  by  damage  entity  ^i.e.,  cohesive  surface)  k.  Shown  in  Fig.  21  is  the  average 
nominal  stress  component  S  =  V~^  fyS^^  dV  for  simulations  2-7,  a  variable  which 
relaxes  towards  zero  magnitude  in  each  case  as  the  applied  deformation  exceeds 
a  =  0.15  and  microcracking  dominates  the  net  elongation  of  the  volume  element. 
Comparing  Figs.  20  and  21,  we  notice  an  inverse  correlation  between  the  volume- 
averaged  nominal  axial  stress  carried  by  the  material,  S  ,  and  damage 
accommodation  in  the  axial  direction,  A\  generally,  the  greater  the  damage 
accommodation  factor  A,  the  greater  the  average  load  reduction  in  the  aggregate 
as  a  result  of  distributed  microcracking. 

For  direct  comparison  with  the  mode  I-type  damage  observed  in  experiments 
(Weerasooriya,  2003)  and  in  our  numerical  calculations,  the  following  evolution 


-  W-matrix  fracture  only  (2) 

.  W-W  fracture  only  (3) 

- W-W  and  W-matrix  fracture  (4) 

—  —  —  -  W-W  and  W-matrix  fracture  (5) 

- W-W  and  W-matrix  fracture  (6) 

- W-W  and  W-matrix  fracture  (7) 


-22 

Fig.  21.  Average  axial  stress,  S  ,  vs.  applied  axial  stretch,  e. 
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equation  is  suggested  below  for  diagonal  components  of  F**: 

iFX=A(F‘‘^-5‘’y,  Fto  =  0,  {a  =  A),  (50) 

where  the  Macaulay  brackets  are  necessary  to  ensure  that  mode  I  cracks  do  not 
contribute  to  compressive  deformation  of  the  volume  element  (i.e.,  no  interpenetra¬ 
tion  of  matter  or  negative  crack  opening  displacements).  Notice  that  Fig.  20  then 
provides  the  time  history  of  A  consistent  with  Eq.  (50).  For  a  given  imposed  total 
deformation  gradient  F,  the  greater  the  magnitude  of  components  of  F*^,  the  less 
strain  accommodation  by  F,  resulting  in  a  reduction  in  magnitude  of  the  average 
stresses  carried  by  material,  as  is  clear  by  comparing  Figs.  16,  20,  and  21. 

Extending  Eq.  (50)  to  arbitrarily  multi-axial  deformations,  we  propose  the 
following  evolutionary  description  for  the  contribution  of  damage  F^  to  the 
deformation  gradient  F,  assuming  an  undamaged,  undeformed  sample  in  the  initial 
state: 


(P)‘:A=5‘'i!'AF%-y^),  F,=o  =  l,  Fto=0,  (51) 

where  is  a  stress-state  dependent  history  variable  of  rank  four — nonzero  upon 
attainment  of  a  critical  nucleation  criterion — scaling  the  degree  of  deformation 
accommodation  from  damage.  Considering  the  effects  of  superposed  rigid  body 
motion  in  the  spatial  frame,  +  c^,  with  the  unimodular  matrix  satisfying 

and  the  spatially-constant  translation  vector  c^,  we  arrive  at  the  following 
objectivity  requirements: 


Cb^.A, 


(F%^aJF‘‘C,  A' 


a.B 

^.b.A 


<,c.B 

.d.A' 


(52) 


Notice  that  Eq.  (50)  is  a  particular  form  of  Eq.  (52),  where  in  Eq.  (50)  we  have 
uncoupled  components  of  evolution  of  damage  deformation  and  total  deformation 
corresponding  to  differing  coordinate  directions  by  letting  A  ^  ^  For 

complex  imposed  multi-axiaMeformations,  relationships  between  components  of  the 
accommodation  measure  A'^^  and  average  stress  reduction  in  each  spatial  direction 
are  expected,  though  specific  forms  of  such  relationships  may  not  be  trivially 
identifiable. 


5.  Conclusions 

Constitutive  models  for  constituents  of  a  two-phase  W-Ni-Fe  tungsten  heavy 
alloy  have  been  developed  within  the  framework  of  continuum  crystal  plasticity 
theory.  The  material  models  capture  finite  deformation,  strain  rate  dependence, 
thermal  expansion,  thermal  softening,  heat  conduction,  and  thermoelastic  heating. 
Dislocation  densities  are  treated  as  internal  state  variables  that  contribute  to  the 
strain  hardening  via  a  square-root  dependency.  Upon  assuming  a  particular  form  for 
the  residual  free  energy  depending  linearly  upon  the  dislocation  density,  the 
framework  allows  computation  of  the  stored  energy  of  cold  working.  A  cohesive 
zone  approach  has  been  used  to  model  fracture  at  W-W  and/or  W-matrix  interfaces 
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in  dynamic  finite  element  calculations  simulating  the  tensile  deformation  of  a 
poly  crystalline  aggregate  of  the  heavy  alloy. 

Our  numerical  results  highlight  the  following  aspects  of  material  behavior  in  the 
WHA  system,  within  the  context  of  high-rate  (a  =  lO'^/s)  tensile  simulations: 

•  W  grains  tend  to  support  higher  stresses  than  the  matrix  phase. 

•  Plastic  strain,  dislocation  accumulation  associated  with  cumulative  strain  hard¬ 
ening,  and  temperature  rise  due  to  plastic  dissipation  are  generally  much  more 
pronounced  in  the  ductile  matrix  phase  than  in  the  initially  stiffer  W  grains. 

•  For  a  uniform  magnitude  of  cohesive  strength  assigned  to  all  interfaces,  fracture 
tends  to  initiate  at  W-W  grain  boundaries  rather  than  at  W-matrix  boundaries. 

•  Formation  of  macro-cracks  of  significant  length  (i.e.,  traversing  the  volume 
element)  requires  activation  of  fracture  sites  at  both  W-W  and  W-matrix 
boundaries. 

•  When  interfacial  fractures  are  simulated,  cracks  accommodate  much  of  the 
deformation  that  would  otherwise  be  accommodated  by  inelastic  straining  in  the 
matrix  phase,  though  regions  of  localized  plastic  flow  and  intense  temperature  rise 
frequently  emerge  in  the  immediate  vicinity  of  damaged  (i.e.,  cohesive)  zones. 

•  Positive  correlation  exists  between  availability  of  fracture  sites,  net  accommoda¬ 
tion  of  deformation  due  to  damage,  and  average  stress  reduction  as  a  result  of 
microcracking. 

Our  simulations  have  also  shown  that  the  macroscopic  effective  stress-strain 
behavior  of  the  aggregate  is  only  mildly  influenced  by  the  grain  arrangement  or  the 
assignment  of  different  random  initial  lattice  orientations  for  either  phase.  However, 
it  remains  to  be  seen  in  future  studies  if  prescription  of  a  preferred  set  of  grain 
orientations  exerts  a  significant  effect  on  the  effective  stiffness  and  the  fracture 
behavior.  Note  that  our  conclusions  are  drawn  in  the  context  of  2D  simulations;  a 
fully  three-dimensional  finite  element  model  of  the  composite  polycrystalline 
microstructure,  with  cohesive  fracture  surfaces  capable  of  tracking  3D  crack 
propagation  and  interaction  (Ortiz  and  Pandolfi,  1999),  would  be  expected  to 
represent  some  improvement  in  terms  of  capability  of  reproducing  experimentally- 
observed  fracture  patterns. 
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